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Abstract 

We formulate a nonlinear continuum model of a graphene sheet supported by a 
flat rigid substrate. The sheet is parallel to the substrate and loaded on a pair of 
opposite edges. A typical cross-section of the sheet is modeled as an elastica. We use 
elementary techniques from bifurcation theory to investigate how the buckling of the 
sheet depends on the boundary conditions, the composition of the substrate, and the 
length of the sheet. We also present numerical results that illustrate snap-buckling 
of the sheet. 
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1 Introduction 

A graphene sheet is a planar hexagonal lattice of carbon atoms with each atom bonded 
to its three nearest neighbors. Theoretical predictions about the properties of graphene 
as well as the role graphene plays as the basic structure in other important materials, in 
particular bulk graphite and carbon nanotubes, have driven efforts to produce isolated 
single-layer graphene sheets. Somewhat surprisingly, these efforts succeeded only recently, 
with the discovery of mechanical and chemical methods for isolating individual graphene 
sheets and functionalized sheets from bulk graphite [23, 24, 36]. The last several years 
have brought a tremendous amount of theoretical and experimental research on graphene. 

Much of this recent research explores the novel electronic transport properties of 
graphene. To a lesser extent, the mechanical properties of graphene have also attracted 
attention. Transport properties suggest engineering nanoscale devices that use graphene 
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as basic components like nanoscale resonators, switches, and valves. See, for example, [3]. 
Mechanical properties suggest the use of graphene in composite materials [27, 39] . For the 
latter, understanding the response of individual graphene sheets to applied loads is clearly 
important. Additionally, for designing graphene-based devices, and more generally for ex- 
ploiting the transport properties of graphene, an understanding of the mechanical response 
of graphene may prove essential for at least two reasons. First, to assemble nanoscale 
components, it may be useful to develop techniques for manipulating individual graphene 
sheets and related nanoscale structures, which entails understanding how these sheets re- 
spond to applied loads [18, 29, 47]. Second, research suggests that the electronic transport 
properties of graphene are coupled to its mechanical deformation [4, 7, 12, 13, 21, 33]. 

In this paper we study the mechanical response of a graphene sheet parallel to and 
supported by a flat rigid substrate. The sheet is loaded compressively on a pair of opposite 
edges. See Figure 1 below. The problem we formulate can be loosely motivated by a 
recent experimental paper of Schniepp et al. [35] , in which functionalized graphene sheets 
supported on a substrate of highly oriented pyrolytic graphite (HOPG) are manipulated 
with the tip of an atomic force microscope (AFM). The authors show that the lateral force 
exerted on the edge of the sheet by the AFM tip can slide the sheet across the substrate 
or, more interestingly, can fold the edge of the sheet. See Figures 3 and 4 in [35]. In fact, 
the sheet can be folded and unfolded by the tip multiple times with this repeated folding 
occurring along the same location. In our idealized continuum model, the compressive load 
on the edges of the sheet could describe the lateral load applied by an AFM tip. More 
generally, the geometry of our problem seems fundamental for the study of the mechanics 
of graphene because mechanical exfoliation, a technique for isolating individual graphene 
from bulk graphite, yields graphene sheets supported by a rigid substrate [24]. Also, 
several proposed electronic devices are based on graphene in a geometry similar to that of 
our problem [33]. 

We consider two types of substrates. In one case, the sheet interacts with SiC>2, a 
material commonly used to isolate single-layer graphene by mechanical exfoliation [2]. The 
interaction between the flexible sheet and the substrate is by van der Waals forces, which 
act over a short-range between the individual carbon atoms on the sheet and the atoms 
in the substrate. Several recent papers suggest that the mechanical response of graphene 
supported on SiC>2 is significantly influenced by the interaction between the sheet and the 
substrate [33]. As a second case, we consider a graphene sheet interacting by van der 
Waals forces with a second, rigid graphene sheet. The interaction between two graphene 
sheets should approximate well the interaction between graphene and the HOPG substrate 
in the experiments in [35] (although we note that the specific interaction terms we present 
in the next section are appropriate for describing a substrate supporting pure graphene, 
not functionalized graphene). 

Modeling both the graphene sheet and the substrate as continua, we derive equilib- 
rium equations for the static response of the compressively loaded sheet. We assume that 
the sheet deforms the same in each cross-section, a simplification supported by the atom- 
istic simulations in [17, 45]. A typical cross-section of the sheet is modeled as an elastic 
rod. Our assumptions along with appropriate boundary conditions — discussed in the next 
paragraph — yield a geometrically exact version of the classical problem of a beam on a 
nonlinearly elastic foundation. 

Two sets of boundary conditions are considered. For the first set, the horizontal com- 
ponent of the applied load at opposite edges of the sheet is prescribed, the loaded edges 
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are moment free, and the loaded edges are kept at a prescribed distance from the sub- 
strate. These correspond to classical 'hinged' boundary conditions. For the second set of 
boundary conditions, we again prescribe the horizontal component of the applied load at 
opposite edges and we require that the loaded edges are moment free. But the loaded edges 
are not kept a prescribed distance from the substrate. Rather, the vertical component of 
the applied load is prescribed to be zero. This second set of conditions is loosely motivated 
by the AFM tip experiments described above, in which the edges of the sheet folded over 
as the lateral load was applied. 

Using standard ideas from bifurcation theory, we investigate how the mechanical re- 
sponse of the graphene sheet depends on the length of the sheet, on the boundary con- 
ditions, and on the composition of the substrate. In our bifurcation analysis, the trivial 
branch corresponds to a flat sheet parallel to and at the equilibrium spacing from the 
substrate. By linearizing about the trivial branch, we compute the critical load at which 
the sheet buckles. How this critical load depends on length, boundary conditions, and the 
composition of the substrate is summarized in Table 1. We show that the critical load 
converges to a nonzero value as the length of the sheet goes to infinity. To illuminate the 
post-buckling behavior, we numerically continue the branching solutions. Our numerical 
results indicate that the post-buckling behavior includes secondary bifurcations, and, more 
interestingly, snap-buckling. 

There is an extensive literature on the buckling of single-walled and multi-walled car- 
bon nanotubes under axial and radial loads. For results that use continuum modeling, 
see, for example, [30, 31, 37]. The literature on buckling of graphene sheets under com- 
pressive loading is less extensive. In [45], we studied the buckling under edge loads of 
two parallel, deformable graphene sheets interacting by van der Waals forces. We showed 
that a simple continuum model of buckling yields predictions that are qualitatively similar 
to the predictions of atomistic simulations. In [17], the authors studied the buckling of 
graphene nanoribbons under compressive load. The atomistic simulations they performed 
were based on energy minimization techniques in which the carbon-carbon interactions 
were modeled by the REBO potential. Their results suggest that linear beam theory pre- 
dicts reasonably well the buckling of graphene. In their problem, the sheet did not interact 
with a substrate. In [34], the author used the 'molecular structural mechanics' approach 
developed in [15] to perform atomistic simulations of the buckling of a graphene sheet 
under compressive loading. How the boundary conditions, the length of the sheet, and 
the aspect ratio of the sheet influenced the buckling load were studied. In their model, 
the sheet is freely suspended and not interacting with a substrate. A comparison be- 
tween Figures 2 and 3 and between Figures 4 and 5 in [34] suggests that for sheets with 
lateral dimensions larger than 25 nm, chirality does not strongly influence the buckling 
load. That the buckling load does not appear to depend on chirality was also noted in 
[17]. In our continuum model, we cannot describe the chirality of the sheet. In [26], the 
authors apply nonlocal elasticity theory to study the buckling of graphene sheets under 
biaxial compression. They study a single graphene sheet that is freely suspended and not 
interacting with a substrate. Figure 5 in [26] suggests that as the lateral dimensions of the 
sheet exceed about 25 nm nonlocal effects become less important for determining buckling 
loads. This observation supports the validity of our results, where nonlocal effects are not 
considered. We note also the recent paper [16], in which the authors study the morphology 
of a graphene sheet supported by a patterned substrate. Their results indicate that as a 
certain parameter controlling the pattern on the substrate is varied, the sheet can exhibit 
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a 'snap-through' instability, which appears similar to the snap-buckling we describe below. 

An interesting issue related to the basic morphology of a graphene sheet supported on 
a substrate is the existence of ripples [19]. Several theories have been proposed to explain 
this rippling [4, 6, 9, 11, 42]. Here we do not directly address the question of rippling. 

In the next section, we derive the boundary-value-problem for a graphene sheet inter- 
acting with a substrate. In Section 3, we identify a trivial branch and compute buckling 
loads. Section 4 contains numerical results that describe post-buckling behavior. The 
final section summarizes our results and mentions several additional problems that could 
be studied within the framework of the model we develop. 

2 Equilibrium Equations 

To describe the basic geometry of our problem, we let {i,j,k} denote a right-handed 
orthonormal basis for IR 3 . The rigid substrate is parallel to the ik-plane, and j points away 
from the substrate. See Figure 1. We assume the deformation of the sheet is the same 
in any cross-section defined by a plane perpendicular to k, and hence the configuration of 
the sheet is determined by the configuration of a typical cross-section. A cross-section is 
described by a curve [0, L] 3s4 r(s) in the ij-plane. 




Figure 1: Basic geometry of a typical cross-section. 

For the geometry depicted in Figure 1, an appropriate theory of nonlinear rods (see 
[1]) delivers the governing equations 

(2.1) n' + f = 0, M' + k- (r' x n) = 0, 

(2.2) n(0) ■ i = n(L) • i = -v, M(0) = M(L) = 0, 

which are the linear and angular momentum balances for the sheet supplemented by bound- 
ary conditions. In (2.1), (2.2), n(s) is the contact force per unit width on the material 
point s, f(s) is the force per unit area exerted by the substrate on s, and M(s) is the k 
component of the contact torque per unit width on s. The boundary condition (2.2) 1 de- 
scribes the load applied to opposite edges of the sheet. The parameter v is the component 
of the load parallel to the substrate at the left edge; v > corresponds to compressive 
loading. The boundary condition (2.2) 2 states that the loaded edges are moment free. 
Additional boundary conditions are given below. 
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To introduce components, we write 



(2.3) 



r(s) = x(s)i+ (y(s) +d)j, 



n(s)=P(s)i + g(s)j, 



where the constant d is defined later. We assume the cross-section of the layer is inexten- 
sible, so that 



for some function [0, L] 9s4 0(s). Hence L is the length of the sheet. We also assume 
that M(s) = a 6'(s) for some constant a > 0, which is a material parameter describing the 
resistance of the sheet to bending. Values from the literature on the continuum modeling 
of graphene and carbon nanotubes suggest that one can choose a in a range between .13 
and .2 nN nm [12, 14, 22, 28, 33, 41, 44, 46]. To make explicit computations below, we 
take a = 0.16 nN nm. Note that inextensibility and a linear dependence of the moment 
on the curvature are the standard assumptions of the elastica theory. 

The term f in (2.1) 1 describes the interaction force exerted by the substrate on the 
sheet. The specific form of / depends on the composition of the substrate. One case we 
study is a substrate composed of SiC>2. To determine the interaction force in this case, we 
start with an expression for the attractive energy between graphene and Si0 2 from [33]. 
The authors report an attractive energy per unit area of 

(,5) *=-£^. 

where g\ = 5.4 • 10~ 3 , g2 = 3.5 ■ 10~ 2 , vf = 10 6 m/s, a = .142 nm , h is Planck's constant, 
and £ is the distance between the graphene sheet and the substrate. From (2.5), we derive 
an attractive force of the form — c^/£ 3 , where ca = 1.499 ■ 10~ 2 nN nm. Now we define 
f (s) = fi(y(s) + dr)j, where // has the form 



with cr = 1.499 • 10~ 2 nN nm 9 . The repulsive term in (2.6) is derived by first assuming a 
dependence, which is consistent with the repulsive part of the Lennard- Jones potential 
used to describe the interaction between non-bonded atoms (see (2.7) below), and then 
choosing cr so that the equilibrium spacing between the sheet and the substrate is dj = 
1 nm, which is consistent with some experimental data [8, 40]. The graph of // is depicted 
in Figure 2(a). 

We also consider the case in which the substrate is HOPG. The force described by f in 
(2.1) 1 arises because of van der Waals interactions between the carbon atoms of the sheet 
and the carbon atoms of the substrate. We consider the interactions between the sheet and 
only the top layer of the atoms on the HOPG substrate. Hence we model the substrate as 
a second, rigid sheet of graphene. (Because the van der Waals force between carbon atoms 
decays rapidly as the spacing between the atoms increases, including interactions with 
additional layers of the HOPG substrate does not significantly change the interaction force 
between the sheet and the substrate.) To define f appropriately for a continuum model, 
we assume the atoms are distributed on the substrate and on the sheet with a uniform 
atomic density a = 38.177 nm~ 2 , a value computed from the geometry of graphene. Also, 
we assume the substrate is infinite in extent. By computing an appropriate improper 



(2.4) 



r'(s) = x'(s)i + y'(s)j = cos6(s)i + sin#(s)j 



(2.6) 



fi(0 
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integral, we find that the force per unit area exerted by the substrate on the sheet is 
f(s) = fn(y(s) + dn)j, where 

(2-7) HO = (p " p) 

with cia = 3.859 x 10" 9 nN nm 13 and c 6 = 2.43 x nN nm 7 . (See [45] for details. The 
values of c 6 and Ci 2 are from [10].) The graph of fjj is depicted in Figure 2(b). Note that 
fn has a unique zero at du = (c^/cq) 1 ^ 6 = 0.341 nm, which is the equilibrium spacing. 

Below we denote the j component of the interaction force by just / and indicate whether 
/ = ji or / = /// only when presenting numerical results. We treat d likewise. Recalling 
(2.3) 1; we see that the j component of r is measured from d, so that y(s) = implies that 
f(a) = 0. 

We now reformulate (2.1), (2.2) in components. First we note that f i = 0, (2.1) l5 and 

(2.3) 2 imply that P'(s) = for all s, and hence P = -v by (2.2) x . Then (2.1), (2.3), and 

(2.4) yield the system 



(2.8a) y' = sin9, 

(2.8b) M' = -z/sin6» -Qcos9, 

(2.8c) & = oT x M, 

(2.8d) o , = -f(v+a). 

We study this system with the boundary conditions 
(2.8e) M(0) = M(L) = 0, 

which is (2.2) 2 , supplemented with either 

(2.8f) y (p) = y(L) = or Q(0) = Q(L) = 0. 
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For conditions (2.8e), (2.8f) 1; the loaded edges are moment free and kept at a prescribed 
distance from the substrate. These conditions correspond to standard hinged boundary 
conditions. For conditions (2.8e), (2.8f) 2 , the edge is moment free but not kept a pre- 
scribed distance from the substrate. Rather, the vertical component of the applied load is 
prescribed to be zero. Below we shall refer to this second set of conditions as the 'floating- 
edge' boundary conditions. (For boundary conditions (2.8f) 1 , the applied load at the edge 
may have a non-zero vertical component.) 
To non-dimensionalize, we define 



(2.9) S=f y(S) = ^, M(S) = ^, 0(g) = 9(s), 

L L a 

(2.10) Q(S) = —Q(s), v=-v, f(0 = -m + d). 

a a a 

Upon inserting these rescalings into (2.8) (and dropping the hats), we get 

(2.11a) y' = sm9, 

(2.11b) M' = —usm9 — Q cos 9, 

(2.11c) 9' = M, 

(2.11d) Q> = -f(y), 

with boundary conditions 

(2.12a) M(0) = M{1) = 0, and 

(2.12b) y ( )=y(l) = or Q(0) = Q(l) = 0. 



Note that the problem we have formulated is a version of the classical problem of a 
beam on an elastic foundation [43]. Our version differs from most treatments because we 
do not use the Euler-Bernoulli beam equation. Also, our choice for /, which describes the 
interaction of the beam with the foundation, is specific to graphene supported by either 
an Si0 2 or HOPG substrate. 

3 Buckling Loads 

One checks that (y, M, 9, Q) = (0, 0, 0, 0) is a solution to (2.11), (2.12) for all loads v. This 
solution corresponds to a flat sheet parallel to the substrate at the equilibrium spacing d. 
In this section, we find the loads at which the sheet buckles from the flat configuration 
and we describe how these buckling loads depend on various parameters in the problem. 
Linearizing (2.11) about the trivial branch yields 



(3.1a) y' = 9, 

(3.1b) M' = -v9 - Q, 
(3.1c) & = M, 

(3.1d) Q' = -f iy , 



where /i := /'(0)_= (L 4 /a)f (d)_by (2.10) 3 . From (2.6), (2.7), we have that /'(d) = -.12 
nN/nm 3 for f — fi and that f'(d) = —84.1 nN/nm 3 for / = ///. The linearized boundary 
conditions for (3.1) are the same as (2.12). 
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To locate critical loads, it is convenient though not essential to reformulate (3.1) as a 
single fourth-order equation. We do so by noting that 



(3.2) 9"" = M'" = -i/0" - Q" = -i/0" + hv = ~v0" + fiO, 
and hence 

(3.3) 9"" + i/0" - j\9 = 0. 

By (3.1c), the condition (2.12a) corresponds to 0'(O) = 0'(1) = 0. Also, by (3.1b) and 
(3.1c), Q = -v0- 9", and then by (3.1d) f\y = v& + 9"'. Hence (2.12) implies that 

(3.4a) 0'(O) = 0'(1) = 0, and 

(3.4b) 0"'(O) = /// (l) = or 0"(O) + 1/0(0) = 0"(1) + i/0(l) = 

are the two sets of boundary conditions for (3.3). 

The critical loads are the values of v at which (3.3), (3.4) has a non-trivial solution. 
For a given u, (3.3), (3.4) has a non-trivial solution if and only if (3.1), (2.12) has a non- 
trivial solution. To find non-trivial solutions of (3.3), (3.4), we compute the characteristic 



roots of (3.3), which are ±y (— v ± ^/ v 2 + 4f\)/2. Using these roots, we find the general 
solution to (3.3), to which we apply one set of boundary conditions from (3.4). This yields 
a system of linear equations M{y)c = for the vector c 6 R 4 of arbitrary constants in the 
general solution. Here M(y) is a 4 x 4 matrix whose entries depend nonlinearly on v. The 
buckling loads are the values of v such that M{y)c = has non-trivial solutions. Hence 
it is sufficient to consider v >->■ det(M(z/)). Performing the computations just described 
yields 

(3.5) det(M(z/)) = sin \ x sin A 2 

for the boundary conditions (3.4a), (3.4b) x and 

(A 6 + A 6 
2(cos A x cos A 2 - 1) + . 1 - . ' sin Ax sin A 2 
{-hr 2 

for the boundary conditions (3.4a), (3.4b) 2 , where 

(3.7) Ai = yV" A 2 = 

The expressions (3.5) and (3.6) are correct for v > 2y/—j\. One can show that there are 
no non-trivial solutions for < v < 2^—f\. Note that f\ and hence det(M(z/)) depend 
on a and L. Also, note that (3.5) implies that either Ai = mix or A 2 = ran, which is 
equivalent to 

(3.8) v = Tx 2 m 2 - fi/n 2 m 2 , 

a formula that could be found more directly by substituting 9 = Acos(rmis) into (3.3), 
where A is an arbitrary amplitude. 

For a given length L, we let v* denote the buckling load, i.e., the smallest critical 
load. In the rest of this section, we illustrate how v* depends on the interaction force, 
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the boundary conditions, the length of the sheet L, and the bending stiffness a. These 
results are computed using either (3.5) or (3.6). For all the results described below, we 
take a = .16 nN nm unless otherwise noted. 

Figure 3 indicates how u* depends on the length L of the sheet for graphene supported 
by an SiC>2 substrate. Figure 3 is for hinged boundary conditions (2.12a), (2.12b)!. From 
the figure we see that as L increases the buckling load oscillates. The local minima of these 
oscillations equal .277 nN/nm, the local maxima decrease, and v* converges to the value 
of the local minima as L gets large. See also Figure 9 (a) below. This oscillatory behavior 
is well-known. See [38] . It occurs because the deformation mode of the buckled solutions 
changes at the cusps in Figures 3 (a), (b). In particular, across the cusps, the nodal 
behavior and the number of half-sine waves in the deformation mode change (see Figure 4, 
which is described in the next paragraph). Each deformation mode has a different critical 
load, and which mode has the lowest critical load changes at the cusps. 




0.276 1 1 1 ' ' 1 ' 

15 20 25 30 35 40 45 50 

L (nm) 
(b) 



Figure 3: Buckling loads for SiC>2 substrate and hinged boundary conditions. Figures (a) 
and (b) show how the buckling load, denoted by u*, depends on L for a = .16 nN nm. 

Figure 4 illustrates how the nodal behavior of the buckled solutions changes at the 
values of L that correspond to the cusps in Figures 3 (a), (b). For the curve describing 
v* as a function of L, there are cusps at L — 4.8, 8.3, 11.7, and 15. The second and third 
of these cusps can be seen in Figure 4 (a). Figure 4 (b) shows the buckled configuration 
corresponding to the buckling load v* for L = 6 nm. The curve in Figure 4 (b) has a single 
simple zero, or node, on (0, L), which is the case for all buckled solutions corresponding 
to lengths L between 4.8 and 8.3. Figures 4 (c) and (d) show how the number of nodes of 
the buckled solutions increases for solutions corresponding to L between 8.3 and 11.7 and 
for solutions corresponding to L between 11.7 and 15. 



9 




Figure 4: Shape of buckled solutions. Figure (a) shows how the buckling load v* depends 
on L for a = .16 nN nm and for graphene supported on an SiC>2 substrate with hinged 
boundary conditions. The points marked I, II, and III correspond to the buckling loads 
for L = 6, 9, and 12 nm. Figures (b), (c), and (d) show how the nodal behavior of the 
corresponding buckled solutions changes across the two cusps located between the points 
labeled I, II, and III in (a). For each solution in (b), (c), and (d), the amplitude has been 
normalized to 1. 



Figure 5 indicates how v* depends on L for graphene supported by an SiC>2 substrate 
and loaded by the floating-edge boundary conditions (2.12a), (2.12b) 2 . As in Figure 3, we 
see that as L increases the buckling load oscillates. In this case, however, both the local 
minima and local maxima decrease as L increases. Figure 5 (a) suggests that v* converges 
to a limiting value as L gets large. See also Figure 9 (a) below. The behavior of the 
buckled solutions across the cusps is similar to that illustrated in Figure 4. 

Figures 6 and 7 illustrate how u* depends on L for graphene supported by a HOPG 
substrate. Figure 6 is for hinged boundary conditions and Figure 7 is for floating-edge 
boundary conditions. Although both Figures 6 and Figure 7 are qualitatively similar to 
Figures 3 and Figure 5, which are the corresponding figures for graphene on SiC>2, we 
note that the predicted buckling loads the HOPG substrate are approximately an order 
of magnitude larger than the buckling loads for graphene supported by an SiC>2 substrate. 
For example, in Figure 3 (b), v* = 0.2771 nN/nm for L = 40 nm while in Figure 6 
(a), v* = 7.335 nN/nm for L = 40 nm. See also Table 1. The difference in the size of 
the buckling loads predicted by the model is a consequence of the difference between the 
derivatives of the interaction forces at the equilibrium spacings — f'i(di) = —.12 nN/nm 3 
for SiO"2 and f' n (du) = —84.1 nN/nm 3 for HOPG. Figure 8 illustrates how u* depends on 
f'(d) for L = 5 nm. Note that —f'(d) represents the linear stiffness or 'spring constant' of 
the elastic foundation. 

The last figure in this section, Figure 9, shows how the limiting value of u* for large L 
depends on a. Figure 9 (a) is for graphene supported by Si02 and for L = 150 nm. As one 
would expect, the figure indicates that for large L the buckling load v* is the same for both 
hinged and floating-edge boundary conditions. The figure also shows that the buckling 
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14 16 18 20 22 24 26 28 30 32 

L (nm) 
(b) 

Figure 5: Buckling loads for graphene supported by an SiC>2 substrate with floating-edge 
boundary conditions. Figures (a) and (b) show how the buckling load v* depends on L 
for a = .16 nN nm. 

load increases as a increases, i.e., as the sheet gets stiffer. Figures 9 (b) corresponds to 
Figures 9 (a) but for graphene supported by a HOPG substrate. 

The results of this section are summarized in Table 1 

4 Post-Buckling Behavior 

In this section we describe the post-buckling behavior predicted by our model. We present 
numerical results computed using the bifurcation software AUTO [5]. 

Figure 10 depicts a part of the bifurcation diagram for a sheet of length L = 30.4 nm 
supported by a HOPG substrate with hinged boundary conditions. (Note that the dashed 
lines are used for clarity and not to indicate that the branches are unstable.) The variable 
on the vertical axis, which measures the size of solutions, is the angle 8(0) between the 
i-axis and the tangent line at the left end of the sheet. The diagram shows a portion of 
the trivial branch containing the buckling load u* and the next two critical loads. There 
is a pitchfork bifurcation at v* . The next two critical loads correspond to transcritical 
bifurcations. We also observe secondary bifurcations on the upper and lower branches 
emanating from the third critical load. Figure 10 should be compared to Figure 11(b) in 
[17] to illustrate the effect of the supporting substrate on the mechanical response of the 
graphene sheet. 

In Figure 11, we illustrate a possible post-buckling path based upon the bifurcation 
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35 40 

L (nm) 
(a) 




(b) 



Figure 6: Buckling loads for HOPG substrate with hinged boundary conditions. Figures 
(a) and (b) show how the smallest eigenvalue, denoted by u*, depends on L for a = .16 
nN nm. 



diagram in Figure 10. We assume the sheet undergoes a quasi-static loading process in 
which a compressive edge load is slowly increased from the unloaded configuration. For 
loads below the buckling load, the trivial branch is stable and the sheet remains flat. We 
assume that at the buckling load the trivial branch loses stability. Because the pitchfork 
opens to the left, if the load is further increased the solution must jump to a different 
branch. For macroscopic structures, this type of jump is often referred to as snap-buckling. 

Determining rigorously to which branch the solution jumps entails analyzing the sta- 
bility of branches by considering an appropriate potential energy for (2.1), (2.2) or by 
studying the dynamical version of these governing equations. Here we do not undertake 
this analysis. (See [32] for a stability analysis based on minimizing potential energy in a 
related problem.) However, to illustrate one possibility, we show the solution jumping to 
the lower branch of the upper pitchfork in Figure 10. The solution follows this branch un- 
til the secondary bifurcation point on the upper branch emanating from the third critical 
load in Figure 10. If the load is further increased, another snap-buckling may occur. The 
shapes of buckled solutions are depicted on the right in Figure 11. Note that both x and 
y are rescaled to be dimensionless. 
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Figure 7: Buckling loads for HOPG substrate with floating-edge boundary conditions. 
Figures (a) and (b) show how the smallest eigenvalue, denoted by u*, depends on L for 
a = .16 nN nm. 

Figure 12 depicts a part of the bifurcation diagram for a sheet of length L = 30.8 
nm supported by a HOPG substrate with hinged boundary conditions. (Here also the 
dashed lines are used for clarity and not to indicate that the branches are unstable.) The 
diagram shows a portion of the trivial branch containing the buckling load v* and the 
next two critical loads. The bifurcation at v* is transcritical. The next two critical loads 
correspond to pitchfork bifurcations. There are secondary bifurcations on the upper and 
lower branches emanating from the second critical load. In Figure 13, we illustrate a 
possible post-buckling path based upon the bifurcation diagram in Figure 12. For loads 
below the buckling load, the trivial branch is stable and the sheet remains flat. At the 
buckling load the trivial branch loses stability. The solution follows the lower branch from 
the transcritical bifurcation until a turning point is reached. If the load is further increased 
the sheet undergoes a snap-buckling and the solution jumps to a different branch. For 
example, the solution could jump to the lower branch of the upper pitchfork in Figure 12. 
The solution follows this branch until the secondary bifurcation point on the upper branch 
emanating from the second critical load. If the load is further increased, another snap- 
buckling may occur. The shapes of buckled solutions are depicted on the right. Both x 
and y are rescaled to be dimensionless. 
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HOPG 



Hinged Floating Hinged Floating 

L min max min max min max min max 

Near 5 nm 0.2770 0.3462 0.4010 0.4617 7.3354 7.4008 7.5417 7.5682 

Near 25 nm 0.2770 0.2795 0.2848 0.2858 7.3354 7.3379 7.3448 7.3456 

Near 150 nm 0.2770 0.2771 0.2774 0.2774 7.3354 7.3355 7.3358 7.3358 



Table 1: Summary of Buckling Loads. The min (max) value corresponds to the buckling 
load at the local minimum (maximum) nearest to the length L on the graph of v* as a 
function of L. The numbers in the table have units of nN/nm. 

The main point of the numerical results in this section is to illustrate the possibility of 
snap-buckling. This possibility suggests that the sheet, if loaded beyond the buckling load, 
would undergo a rapid and relatively large deformation from the flat configuration. As 
noted in the introduction, recent experimental and theoretical work on graphene establishes 
a connection between the deformation of the sheet and its transport properties. Hence our 
modeling predicts that at the buckling load the sheet could undergo a large, rapid change 
in transport properties. This property could be exploited in the design of nanoscale devices 
that incorporate graphene sheets. 

5 Conclusion 

In this paper we developed a nonlinear continuum model of a graphene sheet supported by 
a flat rigid substrate and loaded on a pair of opposite edges. We modeled the cross-section 
of the sheet as an elastica. Using techniques from bifurcation theory, we investigated how 
the buckling of the sheet depends on the boundary conditions, the composition of the 
substrate, and the length of the sheet. We also presented numerical results that illustrate 
some of the possible post-buckling behavior of the sheet. An interesting feature of the post- 
buckling behavior is that the sheet undergoes snap-buckling, which may have implications 
for the design of nanoscale devices that use graphene. 

In our model, we assume the substrate is perfectly flat, which of course is not the 
case. For example, a substrate like Si0 2 may have undulations, and several recent papers 
suggest that graphene supported by an Si0 2 substrate will form ripples to follow these 
undulations [4, 11]. Within the context of our continuum modeling, ripples that form 
in the sheet prior to loading can be described as initial imperfections. A well-developed 
branch of bifurcation theory addresses how the presence of such imperfections influences 
buckling and post-buckling behavior. The study of this question follows naturally from 
the work presented here. 

In a second variation of the problem modeled above, we can assume that the graphene 
sheet is deposited on a deformable, rather than rigid, substrate. Hence the cross-section of 
the substrate could be modeled as a beam with mechanical properties different from the 
sheet. This problem is motivated by recent experimental work in which graphene sheets 
are deposited on substrates of various compositions and strain is induced in the graphene 
by deforming the substrate [20, 21, 25]. 

This work was supported by NSF under grant number DMS-0407361. 
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Figure 8: Buckling loads as a function of the derivative of the interaction force / at the 
equilibrium spacing d. Both hinged and floating-edge boundary conditions are shown. For 
each case, L = 5 nm and a = .16 nN nm. The spacing on the horizontal axis is logarithmic. 
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Figure 9: Buckling loads as a function of a. Figure (a) is for graphene supported by SiC>2 
and Figure (b) is for graphene supported by HOPG. In both cases, L = 150 nm. 



18 



5- o.oon 

-0.005 





* m - m, m 


L = 30.4 nm 






*»», a = .16 nNnm 








) . . / A 

Trivial Branch 












* 






* 



7.330 7J335 7\340 7.345 

v (nN/nm) 



Figure 10: Bifurcation diagram for graphene supported by a HOPG substrate with hinged 
boundary conditions. L = 30.4 nm. The first bifurcation is a pitchfork. The dashed lines 
are used for clarity and not to indicate unstable branches. 
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Figure 11: Possible buckling path. On the left is a part of the bifurcation diagram shown in 
Figure 10. One possible post-buckling path is indicated. On the right are the shapes of the 
buckled solutions along this path. The variables x and y are rescaled to be dimensionless. 
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Figure 12: Bifurcation diagram for graphene supported by a HOPG substrate with hinged 
boundary conditions. L = 30.8 nm. The first bifurcation is transcritical. The dashed lines 
are used for clarity and not to indicate unstable branches. 
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Figure 13: Possible post-buckling path. On the left is a part of the bifurcation diagram 
shown in Figure 12. One possible post-buckling path is indicated. On the right are the 
shapes of the buckled solutions along this path. The variables x and y are rescaled to be 
dimensionless. 
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